Locating the inner edge of neutron star crust using terrestrial nuclear laboratory data 
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Within both dynamical and thermodynamical approaches using the equation of state for neutron- 
rich nuclear matter constrained by the recent isospin diffusion data from heavy-ion reactions in the 
same sub-saturation density range as the neutron star crust, the density and pressure at the inner 
edge separating the liquid core from the solid crust of neutron stars are determined to be 0.040 fm -3 
< p t < 0.065 fm~ 3 and 0.01 MeV/fm 3 < P t < 0.26 MeV/fm 3 , respectively. These together with 
the observed minimum crustal fraction of the total moment of inertia allow us to set a new limit for 
the radius of the Vela pulsar significantly different from the previous estimate. It is further shown 
that the widely used parabolic approximation to the equation of state of asymmetric nuclear matter 
leads systematically to significantly higher core-crust transition densities and pressures, especially 
with stiffer symmetry energy functionals. 
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I. INTRODUCTION 

Having been the major testing grounds of our knowl- 
edge on the nature of matter under extreme condi- 
tions, neutron stars are among the most mysterious ob- 
jects in the universe. To understand their structures 
and properties has long been a very challenging task 
for both the astrophysics and the nuclear physics com- 
munity jl|. Theoretically, neutron stars are expected 
to have a solid inner crust surrounding a liquid core. 
Knowledge on properties of the crust plays an impor- 
tant role in understanding many astrophysical observa- 
tions @, i, 1, i, 1, 0, S i, M EL 03 ■ The inner crust 
spans the region from the neutron drip-out point to the 
inner edge separating the solid crust from the homoge- 
neous liquid core. While the neutron drip-out density 
p ou t is relatively well determined to be about 4 x 10 11 
g/cm 3 [1J], the transition density pt at the inner edge is 
still largely uncertain mainly because of our very limited 
knowledge on the equation of state (EOS), especially the 
density dependence of the symmetry energy, of neutron- 
rich nucleonic matter [y, [7j. These uncertainties have 
hampered our accurate understanding of many impor- 
tant properties of neutron stars [ll H 0] • 

Recently, significant progress has been made in con- 
straining the EOS of neutron-rich nuclear matter us- 
ing terrestrial laboratory experiments (See Ref. [l4[ for 
the most recent review) . In particular, the analysis of 
isospin-diffusion data [la, LLs [l7| in heavy-ion collisions 
has constrained tightly the density dependence of the 
symmetry energy in exactly the same sub-saturation den- 
sity region around the expected inner edge of neutron star 
crust. Moreover, the obtained constraint on the symme- 
try energy was found to agree with isoscaling analyses in 
heavy- ion collisions [18| , the isotopic dependence of the 
giant monopole resonance in even- A Sn isoto pes [19( , and 
the neutron-skin thickness of 208 Pb [IE US IIH . In this 
paper, using the equation of state for neutron-rich nu- 
clear matter constrained by the recent isospin diffusion 



data from heavy-ion reactions in the same sub-saturation 
density range as the neutron star crust, we determine the 
inner edge of neutron star crusts. Consequently, the limit 
on the radius of the Vela pulsar is significantly different 
from the previous estimate. In addition, we find that the 
widely used parabolic approximation (PA) to the EOS 
of asymmetric nuclear matter enhances significantly the 
transition densities and pressures, especially with stiffer 
symmetry energy functionals. 



II. THE THEORETICAL METHOD 

The inner edge corresponds to the phase transition 
from the homogeneous matter at high densities to the 
inhomogeneous matter at low densities. In principle, 
the inner edge can be located by comparing in detail 
relevant properties of the nonuniform solid crust and 
the uniform liquid core mainly consisting of neutrons, 
protons and electrons (npe matter). However, this is 
practically very difficult since the inner crust may con- 
tain nuclei having very complicated geometries, usually 
known as the 'nuclear pasta' 0, OH HI HI, Hi- Fur- 
thermore, the core-crust transition is thought to be a 
very weak first-order phase transition and model calcu- 
lations lead to a very small density discontinuities at the 
transition [f| [25|, [26|, |27J. In practice, therefore, a good 
approximation is to search for the density at which the 
uniform liquid first becomes unstable against small am- 
plitude density fluctuations with clusterization. This ap- 
proximation has been shown to produce very small error 
for the actual core-crust transition density and it would 
yield the exact transition density for a second-order phase 
transition [5J, 1251 . l26l . 1271 ] . Several such methods including 

ip, m m m, the th er - 

3l| and the random phase 
approximation (RPA) [27], l32| have been applied exten 



the dynamical method [2, 
modynamical method 



sively in the literature. Here, we use both the dynamical 
and thermodynamical methods. 



In the dynamical method, the stability condition of 
a homogeneous npe matter against small periodic den- 
sity perturbations with clusterization can be well approx- 
imated by munim 
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ignoring the finite size effects due to surface and Coulomb 
energies as shown in the following. In the above, the 
P = Pb + P e is the total pressure of the npe system with 
the contributions Pi, and P e from baryons and electrons, 
respectively. The v and q c are the volume and charge per 
baryon number. The p is the chemical potential defined 
as 



where 
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and k is the wavevector of the spatially periodic density 
perturbations and pi is the chemical potential of par- 
ticle i. In the above expressions, we used the relation 
dp, „ __ dp T 
dp p ~~ dp„ 

■J^ 2 - with £ being the energy density of the npe mat- 
ter. The three terms in Eq. ((T|) represent, respectively, 
the contributions from the bulk nuclear matter, the den- 
sity gradient (surface) terms and the Coulomb interac- 
tion. For the coefficients of density gradient terms we 
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use their empirical values of D 

MeV-fm 5 consistent with the Skyrme-Hartree-Fock cal- 
culations [H i|. At /fc min = [{^f) 1/2 - k^p] 1 / 2 , the 
Vd yn (k) has the minimal value of Vdy n {k m m) = Vq + 
2(4 7 re 2 /3) 1 / 2 - f3k% F @, i, 0, 1, H|. Its vanishing point 
determines the pt- 

The thermodynamical method re quir es the system to 
obey the intrinsic stability condition [3J] or the following 
inequalities [7l.l30l| 
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In fact, Eq. (|2|) is simply the well-known mechanical sta- 
bility condition of the system at a fixed p. It ensures 
that any local density fluctuation will not diverge. On 
the other hand, Eq. ([3]) is the charge or chemical stability 
condition of the system at a fixed density. It means that 
any local charge variation violating the charge neutrality 
condition will not diverge. If the /3-equilibrium condition 
is satisfied, namely p = p e , the electron contribution to 
the pressure P e is only a function of the chemical poten- 
tial p, and in this case one can rewrite Eq. {2j as 
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By using the relation 
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These conditions are equivalent to require the convexity 
of the energy per particle in the single phase 0, [30| by 



where q c = x p — p e j p. The p = 1/v is the baryon den- 
sity and the E(p, x p ) is the energy per baryon for the 
nucleonic matter. Within the free Fermi gas model, the 
density of electrons p e is uniquely determined by the elec- 
tron chemical potential p e . Then the thermodynamical 
relations Eq. (J2J and Eq. © are identical to @, Hj 
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respectively The second inequality is usually valid. Thus, the following condition from the first one 
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determines the thermodynamical instability region. 
Based on general thermodynamic relations and -g^- = 

-J^, one can show 1331 
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Thus, for d 2 E/dx 2 > 0, Eq. (fT0|) is equivalent to requir- 
ing a positive bulk term Vq in Eq. (|TJ) . Generally speak- 
ing, the pt is in the subsaturation density region where 
d 2 E/dx 2 is almost always positive for all models. There- 
fore, the thermodynamical stability condition is simply 
the limit of the dynamical one as k — > by neglecting 
the Coulomb interaction. 

To locate the inner edge of neutron star crust, we use 
the same MDI (Momentum Dependent Interaction) EOS 
that was used in analyzing the isospin diffusion in heavy- 
ion reactions [lg, LLZ| ■ The MDI interaction is based on a 
modified finite-range Gogny effective interactio n 13511 and 
has been extensively used in our previous work [14J . The 
baryon potential energy density part of the MDI EOS 
can be expressed as [la, bis [36| 
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Here r is 1/2 ( — 1/2) for neutrons (protons) and 6 = 
1 — 2x p is the isospin asymmetry. The meaning and val- 
ues of other parameters can be found in Refs. [la, l3a |. 
The parameter x is introduced to vary only the density 
dependence of the symmetry energy while keeping other 
properties of the nuclear EOS fixed [la]. In particular, 

the symmetry energy E sym (p) = ^ I 4J ) at satura- 

V / 5—0 

tion density po = 0.16 fm~ 3 is fixed at 30.54 MeV for all 
values of the parameter x. The isospin symmetric part of 
the MDI EOS was shown to agree with the experimental 
constraints obtained from relativistic heavy-ion collisions 
up to about 5 times the saturation density [37| ■ 



III. RESULTS 

Firstly, to show how the uniform npe matter becomes 
stable from unstable with increasing baryon density and 
how to locate the transition density and see the difference 
between the dynamical and thermodynamical methods as 
well as effects of the PA, we show in Fig. [T] the density 
dependence of Vd yn and V^ her using the MDI interaction 
with x = within both the dynamical and thermody- 
namical methods with the full EOS and its PA. Here, we 
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FIG. 1: (Color online) The density dependence of Vdyn and 
Vf her for MDI interaction with x = using both the dynami- 
cal and thermodynamical methods with the full EOS and its 
parabolic approximation (PA). 
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and it should be noted that V( h 



has the same vanish- 
ing point as the Vther and the same dimension as the 
Vdyn- For the MDI interaction with x = the transi- 
tion densities using the full EOS within the dynamical 
and thermodynamical method are 0.065 fm~ 3 and 0.073 
fm~ 3 , respectively. While the corresponding results us- 
ing the PA are 0.080 fm~ 3 and 0.090 fm~ 3 , respectively. 
Thus, the transition densities are generally lower with 
the dynamical method as the density gradient term and 
the Coulomb interaction make the system more stable. 
However, the PA significantly lifts the transition density 
regardless of the approach used. In fact, the difference 
between calculations using the full EOS and its PA is 
much larger than that caused by using the two different 
methods. 

Shown in Fig. [5] is the pt as a function of the slope 
parameter of the symmetry energy L = 3po — s g \p=po 
with the MDI interaction. For comparisons, we have in- 
cluded results using both the dynamical and thermody- 
namical methods with the full EOS and its parabolic ap- 
proximation, i.e., E(p 7 S) = E(p,S = 0) + E sym (p)S 2 + 
0(S 4 ) from the same MDI interaction. With the full MDI 
EOS, it is clearly seen that the pt decreases almost lin- 
early with increasing L for both methods. This feature 
is consistent with the RPA results j32[. It is interest- 
ing to see that both the dynamical and thermodynam- 
ical methods give very similar results with the former 
giving slightly smaller p t than the later (the difference 
is actually less than 0.01 fm -3 ) and this is due to the 
fact that the former includes the density gradient and 
Coulomb terms which make the system more stable and 
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FIG. 2: (Color online) The pt as a function of L from the 
dynamical and thermodynamical methods with and without 
the parabolic approximation in the MDI interaction. The 
triangles are obtained by Kubis [30I | and the star with error 
bar represents L — 86 ± 25 MeV. 



lower the transition density. The small difference be- 
tween the two methods implies that the effects of density 
gradient terms and Coulomb term are unimportant in de- 
termining the p t . On the other hand, surprisingly, the PA 
drastically changes the results, especially for stiffer sym- 
metry energies (larger L values) . Also included in Fig. [3] 
are the predictions by Kubis using the PA of the MDI 
EOS in the thermodynamical approach [3fJ. Further- 
more, we find that in the parabolic approximation, using 
the EOS of E(p, 5) = E(p, S = 0) + [E(p, 5 = 1)-E(p, 5 = 
0)](5 2 + 0(5 4 ) leads to almost the same results. The large 
error introduced by the PA is understandable since the (3- 
stablc npe matter is usually highly neutron-rich and the 
contribution from the higher order terms in S is apprecia- 
ble. This is especially the case for the stiffer symmetry 
energy which generally leads to a more neutron-rich npe 
matter at subsaturation densities. In addition, simply 
because of the energy curvatures involved in the stabil- 
ity conditions, the contributions from higher order terms 
in the EOS are multiplied by a larger factor than the 
quadratic term. These features agree with the early find- 
ing [38[ that the p t is very sensitive to the fine details of 
the nuclear EOS. We notice that the EOS of asymmet- 
ric nuclear matter always contains the higher-order terms 
in isospin asymmetry (at least for the kinetic part of the 
EOS). Our results indicate that one may introduce a huge 
error by assuming a priori that the EOS is parabolic for 
a given interaction in calculating the p t . We thus apply 
the experimentally constrained L to the p t — L correla- 
tion obtained using the full EOS in constraining the pt- 
In the following, we will mainly focus on the dynamical 
method as it is more complete and realistic. 

The transport model analysis of the isospin diffusion 
data from heavy-ion collisions allowed us to constrain the 



parameter x in Eq. (|12[) to be between x = and x = — 1 
in the density range of 0.3po an d 1.2p 0, LL21- With 
the full MDI EOS, the slope parameter was determined 
to be L = 86 ± 25 MeV. The approximately 30% error in 
L is systematic in nature mainly due to the uncertainty 
of the in-medium nucleon-nucleon cross sections used in 
the transport model [l7[- The statistical errors in both 
the data analysis [I5[ and model calculations [lj| LLZ[ are 
less than the systematic error. The error in L will lead 
to roughly similar systematic errors in all quantities we 
study here. As shown in Fig. [2j the constrained L then 
limits the transition density to 0.040 fm~ 3 < p t < 0.065 
fm~ 3 . It is interesting to mention that both approaches 
used here for finding the core-crust transition density in 
the npe matter in neutron stars have also been widely 
used in studying the spinodal decomposition density as- 
sociated with the liquid-gas phase transition in nuclear 
matter, see, e.g., Refs. [33, [4fJ, |4l], [42| for some recent 
applications in neutron-rich matter. As expected, see, 
e.g., Ref. [43|, the two phase transitions are asymptoti- 
cally but inherently related. Applying both the dynam- 
ical and thermodynamical approaches to symmetric nu- 
clear matter (SNM) at zero temperature for MDI inter- 
action by ignoring the Coulomb and surface terms, we 
find a transition density of 0.63po- It is consistent with 
the spinodal decomposition density for SNM at zero tem- 
perature from the relativistic mean field model [39( and 
Skyrme density functionals [41j. On the other hand, re- 
cent studies J4J, [45| indicate that some effects beyond the 
mean-field approximation may affect the EOS of asym- 
metric nuclear matter in the low-density region. Thus, it 
would be interesting to study how the effects beyond the 
mean-field approximation as well as other effects such as 
the many-body forces may change the core-crust transi- 
tion density of neutron stars. 
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FIG. 3: (Color online) The P t as functions of L and pt by 
using the dynamical method with and without parabolic ap- 
proximation in the MDI interaction. The star with error bar 
in the left panel represents L — 86 ± 25 MeV. 

The pressure at the inner edge, P tl is also an impor- 
tant quantity which might be measurable indirectly from 
observations of pulsar glitches [3, 0]. Shown in Fig. [3] 
is the Pt as functions of L and pt by using the dynam- 



ical method with both the full MDI EOS and its PA. 
Again, it is seen that the PA leads to huge errors for 
larger (smaller) L (p t ) values. For the full MDI EOS, the 
Pt decreases (increases) with the increasing L (p t ) while 
it displays a complex relation with L or p t for the PA. 
The complex behaviors are due to the fact that the pt 
does not vary monotonically with L for the PA as shown 
in Fig. [5] From the constrained L values, the Pt is limited 
between 0.01 MeV/fm 3 and 0.26 MeV/fm 3 . 

The constrained values of p t and Pt have important 
implications on many properties of neutron stars 0, @, 
lid. l28||. As an example, here we examine their impact on 
constraining the mass-radius (M-R) correlation of neu- 
tron stars. The crustal fraction of the total moment of 
inertia AI/I can be well approximated by p, 0, |9J 
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where rrib is the mass of baryons and £ = GM/Rc 2 with 
G being the gravitational constant. As it was stressed 
in Ref. p, the AI/I depends sensitively on the symme- 
try energy at subsaturation densities through the P t and 
pt, but there is no explicit dependence upon the higher- 
density EOS. So far, the only known limit of AI/I > 
0.014 was extracted from studying the glitches of the Vela 
pulsar |9(. This together with the upper bounds on the 
Pt and p t (p t = 0.065 fm" 3 and P t = 0.26 MeV/fm 3 ) sets 
approximately a minimum radius of R > 4.7 + 4.QM/M & 
km for the Vela pulsar. The radius of the Vela pulsar 
is predicted to exceed 10.5 km should it have a mass of 
IAM@. A more restrictive constraint will be obtained 
from the lower bounds of p t = 0.040 fm~ 3 (P t = 0.01 
MeV/fm 3 ) which is indicated by the curve with solid stars 
in Fig. [4] and it can be approximately parameterized by 
R = 5.5 + 14.5A//A/ km. It is thus seen that the er- 
ror of the transition density and pressure obtained in the 
present work is still large and it leads to large uncer- 
tainties for the mass-radius relation of the Vela pulsar. 
As a conservative estimate, we thus deduce a constraint 
of R > 4.7 + 4.OM/M km using the upper bounds on 
the Pt and pt obtained in the present work. We notice 
that a constraint of R > 3.6 + 3.9M/M km for this 
pulsar has previously been derived in Ref. [9( by using 
p t = 0.075 fnT 3 and P t = 0.65 MeV/fm 3 . However, 
the constraint obtained in the present work using for the 
first time data from both the terrestrial laboratory ex- 
periments and astrophysical observations is significantly 
different and actually it is more stringent. 

To put the above constraints on the Vela pulsar in 
perspective, we show them in Fig. 0] together with 
the M— R relation by solving the Tolman-Oppcnheimer- 
Volkoff (TOV) equation. In the latter, we use the well- 
known BPS EOS [2J for the outer crust. In the inner crust 
with p out < p < pt, the EOS is largely uncertain and fol- 
lowing Ref. 
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FIG. 4: (Color online) The M-R relation of static neutron 
stars from the full EOS and its parabolic approximation in 
the MDI interaction with x — and x = — 1. For the Vela 
pulsar, the constraint of AI/I > 0.014 implies that allowed 
masses and radii lie to the right of the line linked with solid 
squares (p t = 0.065 fm" 3 and P t = 0.26 MeV/fm 3 , the upper 
bounds obtained in the present work), solid stars (pt = 0.040 
fm~ 3 and P t = 0.01 MeV/fm 3 , the lower bounds obtained 
in the present work) or open squares (pt = 0.075 fm -3 and 
Pt = 0.65 MeV/fm 3 , used in Ref. 0). 



with the constants a and b determined by the total pres- 
sure P and total energy density e at p ou t and pt ■ The full 
MDI EOS and its parabolic approximation with x = 
and x = — 1 are used for the uniform liquid core with 
p > pt- In this way, the P is a continuous function of the 
e at the boundary between different regions as required. 
Assuming that the core consists of only the npe matter 
without possible new degrees of freedom or phase tran- 
sitions at high densities, the PA leads to a larger radius 
for a fixed mass compared to the full MDI EOS. Further- 
more, using the full MDI EOS with x — and x = — 1 
constrained by the heavy-ion reaction experiments, the 
radius of a canonical neutron star of 1.4M is tightly 
constrained within 11.9 km to 13.2 km, which is consis- 
tent with the constraint R > 4.7 + 4.0M/M o km for the 
Vela pulsar. 



IV. SUMMARY 

In summary, the density and pressure at the inner edge 
separating the liquid core from the solid crust of neutron 
stars are determined to be 0.040 fm -3 < p t < 0.065 fm -3 
and 0.01 MeV/fm 3 < P t < 0.26 MeV/fm 3 , respectively, 
using the MDI EOS of neutron-rich nuclear matter con- 
strained by the recent isospin diffusion data from heavy- 
ion reactions in the same sub-saturation density range 
as the neutron star crust. These constraints allow us 
to set a new limit on the radius of the Vela pulsar. It 
is significantly different from the previous estimate and 



thus puts a new constraint for the mass-radius relation of 
neutron stars. Furthermore, we find that the widely used 
parabolic approximation to the EOS of asymmetric nu- 
clear matter leads systematically to significantly higher 
core-crust transition densities and pressures, especially 
for the energy density functional with stiffer symmetry 
energies. Our results thus indicate that one may intro- 
duce a huge error by assuming a priori that the EOS is 
parabolic with respect to isospin asymmetry for a given 
interaction in locating the inner edge of neutron star 
crust. 
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